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We investigate the phase structure of the Nambu-Jona-Lasinio model at zero tempera- 
ture, allowing for a two-dimensional spatial dependence of the chiral condensate. Applying 

CNJ ' the mean-field approximation, we consider various periodic structures with rectangular and 

hexagonal geometries, and minimize the corresponding free energy. We find that these two- 

^Sl . dimensional chiral crystals are favored over homogeneous phases in a certain window in the 

region where the phase transition would take place when the analysis was restricted to homo- 
geneous condensates. It turns out, however, that in this regime they are disfavored against 
a phase with a one-dimensional modulation of the chiral condensate. On the other hand, we 

^^ . find that square and hexagonal lattices eventually get favored at higher chemical potentials. 

Although stretching the limits of the model to some extent, this would support predictions 
f~| ' from quarkyonic-matter studies. 
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INTRODUCTION 



Since more than three decades, the phase diagram of quantmn chromodynamics is object of 
intensive theoretical and experimental research prt^]. In particular the conjecture that at low 
temperatures there could be a first-order chiral phase transition which ends at a critical point has 
received a lot of interest, since this endpoint is potentially detectable in heavy-ion experiments. 
However, in most studies which support this scenario, it is tacitly assumed that the order parame- 
ters of the various phases are uniform in space. On the other hand, it is not a new idea that there 
could be spatially modulated states in strongly interacting matter (see Ref . [5| for a recent review) . 
Well-known examples are the proposal of an inhomogeneous ground state in nuclear matter [6|, the 
possibility of spatial modulations in the context ofpion condensation [7|, la], Skyrme crystals [9|], 
and crystalline phases in (color-) superconductors 1CH17I|. 

For the QCD phase diagram, it has been argued some time ago that, at least in the limit of 
a large number of colors (Nc), the favored ground state of a dense Fermi sea of quarks should be 
characterized by a spatial modulation of the chiral condensate [l8|, ll9|. More recent studies on 
quarkyonic matter seem to support this hypothesis 2014221]. 

For the physical case of three colors, Nambu-Jona-Lasinio- (NJL-) type model studies have 
revealed at intermediate chemical potentials and low temperatures the presence of an inhomoge- 
neous phase where the chiral condensate assumes a spatially modulated form [2314251] . Most of the 
existing studies on inhomogeneous phases have restricted their analysis to simplified shapes of the 
chiral order parameter. The most popular example is the so-called "chiral density wave" (CDW), 
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261 ]. This kind of ansatz is analogous to 



which basically amounts to a single plane wave 

the so-called Fulde-Ferrel solutions proposed in (color-) superconductivity [IC 

The study of a generic shape for the spatially dependent chiral condensate is a highly non- 
trivial task. It has recently been observed, however, that in NJL-type models the evaluation of the 
energy spectrum simplifies considerably when the condensate is allowed to vary only in one spatial 
dimension, while remaining constant along the two transverse directions 25|. In this particular case, 
the problem can be reduced to the 1-l-l-dimensional chiral Gross- Neveu model, where analytical 
expressions for the eigenvalue spectra are known 
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301 ]. This formal resemblance allows to perform 



an analysis of the phase diagram including these inhomogeneous phases without having to calculate 
numerically the energy spectrum of the model. 

Within this framework it has been found that the inhomogeneous phase covers the region 
where a first-order chiral phase transition would occur when limiting the analysis to homogeneous 
phases. As a consequence the chiral critical point disappears from the phase diagram, leaving only 
a Lifshitz point where three second-order lines meet 25|, l31[. The inclusion of vector interactions 
further enhances this effect and enlarges the size of the inhomogeneous phase [32] ■ 

The limitation to one-dimensional structures is of course a strong one. Especially at lower 
temperatures, hi ghe r dimensional modulations are expected to play an important role in (color-) 
superconductors J12l.ll4|. whose dynamics bears strong formal resemblance with the one described 
in the NJL model. Moreover, recent quarkyonic-matter studies suggest that in the high chemical 
potential region, as the density of the system grows, the quark Fermi sea tends to break chiral 



symmetry by forming increasingly complex crystalline (or quasi-crystalline) structures, which can 
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be described as superpositions of several "quarkyonic chiral spirals" 

Aside from these considerations, modulations in more than one spatial dimension are also of 
interest since they are unaffected by the instabilities with respect to fluctuations which prevent the 
formation of a true ID crystalline structure at finite temperature 33l |. 

The main purpose of this paper is therefore to investigate the properties of modulations of the 
chiral condensate occurring in more than one spatial dimension. While a complete analysis would 
in principle have to allow for modulations in three spatial dimensions, for computational reasons we 
will limit ourselves in the present work to two-dimensional crystalline shapes. The generalization 
to three-dimensional structures is straightforward once the formalism and the numerical procedure 
are set up. 



II. INHOMOGENEOUS PHASES IN THE NJL MODEL 



Our starting point is the two-flavor NJL Lagrangian 34i |. 



(1) 



where ip \s a. quark field with two flavor and three color degrees of freedom and bare mass m, Ta 
denotes the three Pauli matrices in isospin space, and G is a coupling constant. 

We perform the mean-field approximation by expanding the interaction around the scalar and 
pseudoscalar condensates 



(VV') = S{X) , {i^ij^T'^tP) = P{x) 6a3 , 



(2) 



which we allow to be space dependent. For later convenience, we also introduce the complex "mass" 
function, 



M{x) =m- 2G{S{x) + iP{x)) . 
The mean-field Lagrangian can then be written as 

Cmt = i^i\ido - n)i^ - G{S^ + P2) 

with the effective Hamiltonian operator 



'H = j'' 



i^-d + m- 2G{S + ij^TsP) 



(3) 



(4) 



(5) 



The mean-field thermodynamic potential per volume 0, associated with these (so far generic) spa- 
tially modulated condensates at temperature T and quark chemical potential fi contains a func- 
tional trace over the logarithm of the inverse quark propagator [35| . For its evaluation we employ 
imaginary time formalism and switch to momentum space. Assuming static (i.e., time- independent) 
condensates, we can perform the sum over Matsubara frequencies explicitly and obtain (up to a 
constant) 



0(T, fi- Mix)) = n^niT, fi- M{x)) + ftcondiMix)) 



(6) 



with 



^cond{M{x)) = - j ^'^^ ^> (7) 

where V is the volume of the system, and 

f^,,„(r, ^; M{x)) = -TY, log (2 cosh (^^^) ) , (8) 

where the sum runs over all eigenvalues ii^ of "^^ in color, flavor, Dirac and momentum space. 

In presence of an inhomogeneous condensate, the diagonalization of "^^ is a highly non-trivial 
task, since quarks may exchange momenta by scattering off the condensate and, consequently, the 
resulting mean-field quark propagator is not diagonal in momentum space. In the following we 
assume a periodic shape of the chiral condensate forming a well-defined lattice structure. This 
implies that we can expand the spatially varying order parameter in a Fourier series, 

M(x) = ^M,-,e^^^-^ (9) 

with discrete momenta qj. forming a reciprocal lattice (RL). A generic element oi% in momentum 
space then takes the form 



\ l^qh ^"qk P^,Pn-qk " ■ Pm "pm,Pn / 

where Y^^ runs over the momenta of the RL, making obvious the non-diagonal momentum struc- 
ture of the matrix|^ In turn, momenta which do not differ by an element of the RL are not coupled, 
so that Ti can be decomposed into a block diagonal form, where each block ^{(k) can be labelled 
by an element of the first Brillouin zone (BZ). This implementation of Bloch's theorem allows to 
decompose the eigenvalue sum in Eq. ([8]) into a momentum integration over the BZ times a sum 



over the discrete eigenvalues of each block [17 1 



III. TWO-DIMENSIONAL MODULATIONS 

For general periodic structures, although the numerical diagonalization procedure is in principle 
straightforward, its practical implementation turns out to be computationally demanding. In 
order to simplify the problem, we therefore limit the generality of our ansatz Eq. Q to lower- 
dimensional modulations. In this case, the momentum integration required for the evaluation of 
the thermodynamic potential may be split into the parts pn along the direction of the modulation 
and p_j_ perpendicular to it. One thus obtains, for a d-dimensional modulation, 



d^^'^px r d'^k 



^^- = -y (2^0^/(2^^^^°^ 



BZ ^ 



2co,hp^^ 



(11) 



^ The matrix is also non-diagonal in Dirac space, as indicated in Eq. (|10|) . Here the chiral representation was used, 
and a corresponds to the Pauli matrices. On the other hand, "H is diagonal in isospin and color. 



where k labels the BZ momenta and E{p±, k) are the eigenvalues of T-L{k) for a given p±. 

It has been observed that the eigenvalues in 3+1 dimensions can be evaluated by diagonal- 
izing a dimensionally reduced 7i (evaluated at p± = 0) and subsequently boosting the resulting 
spectrum \23\. This dramatically simplifies the calculations. In particular for the case of one- 
dimensional modulations it allows to reuse a well-established set of analytical results without hav- 
ing to resort to a numerical diagonalization of the model Hamiltonian. The favored mass functions 
are then given by Jacobi elliptic functions, which smoothly interpolate between solitonic shapes 
close to the homogeneous chirally broken phase and sinusoidal shapes close to the restored phase. 

Since the one-dimensional problem has already been treated extensively in previous works, the 
main focus of this work is on two-dimensional structures. A recent Ginzburg-Landau (GL) analysis 
has shown that close to the Lifshitz point one-dimensional modulations are energetically favored 



over higher dimensional ones in NJL-type models [Sg]. We therefore focus on what happens at 
zero temperature, where GL arguments are unable to provide reliable results. For simplicity, we 
restrict our calculations to the chiral limit, m = 0. 

Without loss of generality, we assume the chiral condensate to vary in the xy-plane and to be 
constant in the z-direction. Unlike for the one-dimensional case, in two spatial dimensions different 
crystalline shapes may be realized. Therefore we consider different lattice geometries and assume 
the mass functions to have simple symmetric shapes consistent with these structures. 

The first case is a square lattice with a unit cell spanned by two perpendicular vectors of length a 
in X and y direction. The corresponding elements of the RL are then given by qm,n = Q{mex + ney) 
with Q = 27r/a and integers m and n. While the general mass function consistent with this lattice 
structure would be given by Eq. ([9]) with arbitrary Fourier coefficients Mm,ni we restrict ourselves 
to a simple ansatz for a real symmetric mass function with a small number of nonvanishing Fourier 
coefficients. Specifically we choose Mi^i = Mi__i = M_i^i = M_i^_i = M/4 and Mm,n = in all 
other caseso This yields 

M(x, y) = M cos(Qx) cos{Qy) , (12) 

which has an egg-carton-like shape (see Fig. [H left) and is symmetric under discrete rotations by 
7r/2. 

The second case we consider is a mass function with hexagonal symmetry. Here we start from 
a unit cell spanned by two vectors of length a, enclosing an angle of vr/S. Choosing the first one 
to be aligned with the x-axis, the elements of the RL are given by qm,n = Qifncx + ^"""^ e^y), 
with integers m and n, and Q = 2'ir/a as before. For the mass function we choose Mm,n = M/6 
on the corners of a regular hexagon, {m,n) G {(1,0), ( — 1, 0), (0, 1), (0, — 1), (1, 1), (— 1, — 1)}, and 
Mm,n = in all other cases. This yields 



M{x,y) = — 



2 cos (Qx) cos I —i=Qy ) + cos{—j=Qy) 



(13) 



An alternative choice would be Mi,o = Af_i,o = Mo,i = Mo,_i = M/4 and Mm,n = in all other cases, which 
yields M{x,y) = ^(cos(Q2;) + cos{Qy)). However, this is equivalent to Eq. (|T2)| in a frame rotated by n/4 and Q 
replaced by Q/\/2. 





FIG. 1. Mass functions M{x,y) with two-dimensional modulations in coordinate space. Left: "egg-carton" 
modulation on a square lattice, Eq. (IT2l) . Right: hexagonal modulation, Eq. p3|. 
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FIG. 2. Amplitude M and wave number Q at T = as functions of the chemical potential jj, after minimiza- 
tion of the thermodynamic potential for given shapes of the mass function. Left: "egg-carton" modulation 
on a square lattice, Eq. ([T2l) . Right: hexagonal ansatz, Eq. (|13|) . 



which is symmetric under discrete rotations by 7r/3 (see Fig. [H right). Note that the normalization 
of the amplitude was chosen to match the homogeneous case M(x, y) = M when Q goes to zero. 

After inserting these shapes into Eq. (jlOp . the Hamiltonian is diagonahzed numerically in Dirac 
and momentum space, and the thermodynamic potential is minimized with respect the variational 
parameters M and Q, characterizing amplitude and period of the modulations. For the numerical 
calculations, we have to specify a procedure to regularize the integrals and eigenvalue sum in 



Eq. (jlip . Following Refs. 25l. |32|. |37| . we use a Pauli-Villars scheme with three regulators. For the 



results in this section, we fit our model parameters (the regulator A and the coupling constant G) 
to reproduce the pion decay constant in the chiral limit, /^r = 88 MeV, and a vacuum constituent 
quark mass of M^ac = 300 MeV [zsj. 

The results of the numerical minimization for the square and hexagonal shapes, Eqs. (J12p and 
p^ respectively, are presented in Fig. [2j For both modulations, we find a sharp onset of the 
crystalline phase around // ~ 310 MeV and a smooth approach to the restored phase, which 
is reached at /i ~ 345 MeV via a second-order phase transition as the amplitude of the chiral 
condensate melts to zero. For the transition to the homogeneous chirally broken phase, the situation 
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FIG. 3. Thermodynamic potential relative to the restored phase for different modulations of the chiral 
condensate at T = 0. The homogeneous broken and restored solutions are disfavored compared to the 
crystalline phases in a window between /i « 308 and fi « 345 MeV. The lowest free energy is found for the 
one-dimensional Jacobi elliptic function. In particular at the onset of the inhomogeneous phase it is favored 
over a one-dimensional cosine, while with increasing /i the two shapes quickly become almost degenerate. 
The two-dimensional square ansatz, Eq. ()12p . is always disfavored against the one-dimensional modulations, 
and the hexagonal crystal, Eq. (jl3|) . leads to an even smaller gain in free energy. 



is therefore different from the case of one-dimensional solitonic solutions [25] and might be due to 
our limited ansatz with a finite number of Fourier components. We also note that the results for 
the amplitude M are comparable for both shapes in the inhomogeneous window. 

The results presented in Fig. [2] have been obtained by enforcing in each case a fixed shape of the 
chiral modulation. Under this restriction we found a window where the different two-dimensional 
solutions are energetically favored over the homogeneous solutions. This is basically the same region 
where one-dimensional modulations are also found to be favored over homogeneous solutions. In 
fact, as shown in Ref. 3l[ using GL arguments, the critical chemical potential for the transition 
from the inhomogeneous to the chirally restored phase is independent of the shape of the spatial 
ixLodulation, as long as the phase transition is second order. 

The next obvious step is to compare the free energies of these solutions with each other, in order 
to find out which of them corresponds to the most stable solution. The results of our comparison 
are shown in Fig. [3l One can clearly see that the one-dimensional Jacobi elliptic functions lead 
to the biggest gain in free energy compared to all the other cases considered. In particular, the 
two-dimensional structures turn out to be energetically disfavored with respect to one-dimensional 
real modulations throughout the whole inhomogeneous window. 

In this context it is instructive to introduce a rectangular structure, which interpolates contin- 
uously between a square lattice and a one-dimensional periodic modulation. Specifically, we can 
generalize the "egg-carton" ansatz Eq. (I12p to 

M(x, y)=M cos(Q^x) cos{Qyy) , (14) 

which reduces to a single cosine varying in one spatial dimension when one of the two wave numbers 
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FIG. 4. Value of the thermodynamic potential at T = and ^ = 325 MeV for a rectangular lattice, Eq. p^ . 
as a function of Qj. and Qy. At each point, J7 was minimized with respect to M . 



goes to zero. 

Starting from this ansatz, we minimize the thermodynamic potential with respect to the am- 
plitude M for fixed wave numbers Qx and Qy, and then study the result as a function of Qx and 
Qy. Since we know already that the square-lattice solution along the line Qx = Qy is disfavored 
against the one-dimensional cosine, we are mainly interested in the question whether the former 
corresponds to a local minimum or to a saddle point in the Qx — Qy plane. 

Our result for /i = 325 MeV is presented in Fig. IU showing that the solution at Qx = Qy is 
a local minimum. One can also see that in the vicinity of the minimum the potential remains 
rather flat along the direction perpendicular to the line Qx = Qy Going along this valley, we find 
two saddle points at iQx,Qy) ~ (175 MeV, 400 MeV) and (400 MeV, 175 MeV). Unfortunately, the 
computing time rises strongly with decreasing wave numbers, so that we could not continue this 
analysis to values of Qx or Qy lower than 100 MeV. However, it is not hard to imagine how beyond 
the saddle points the valley approaches the absolute minima at vanishing Qx or Qy, corresponding 
to a one-dimensional cosine. 

One may ask whether our observation that the two-dimensional crystalline structures are dis- 
favored against the one-dimensional ones is caused by the restricted ansatz for the mass functions. 
Taking into account more Fourier modes would lead to additional variational parameters and could 
thus lower the free energy. It is however unlikely that this would change our results considerably. 
As seen in Fig. [3] the difference in free energy between the Jacobi elliptic function and the one- 
dimensional cosine is negligible in a large chemical-potential range and always much smaller than 
the difference to the two-dimensional solutions. We therefore expect that the corrections to the 
considered two-dimensional shapes are small as well. 



In order to test this, we extend the simple "egg carton" ansatz (Eq ll2p into 
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M{x, y) = 2. ^n cos{nQx) cos{nQy) (15) 

n=l 

The minimization of the thermodynamic potential with respect to the variational parameters 
(Mi,M2,M3, Q) leads, within numerical errors, to M2 = M3 = throughout the whole inho- 
mogeneous window. Although numerical errors are more significant than for the one-dimensional 
case, it is safe to say that the inclusion of those higher harmonics considered above does not lead 
to an appreciable gain in free energy. 

IV. HIGHER CHEMICAL POTENTIALS 

In recent quarkyonic-matter studies it was found that increasing the chemical potential leads to 
two-dimensional structures with growing geometrical complexity |22f]. While no definite scale was 
given, this could be a hint that the chemical potentials we have considered so far are too low for 
two-dimensional crystals to be favored. This motivates us to extend our investigations to higher 
values of /x. 

In fact, while so far we have concentrated on the inhomogeneous "island" close to the would- 
be first-order phase boundary for homogeneous phases, we have recently found that in the NJL 
model a second inhomogeneous region ("continent") appears at higher /i and seemingly persists to 
arbitrarily high chemical potentials [37|. Of course, we have to keep in mind that the NJL model 
is a low-energy effective model with a limited range of validity. In particular, since the continent 
appears in a region where the chemical potential is of the order of the regulator masses, we have 



to be cautious not to overinterpret the results. However, as thoroughly discussed in Ref. |37l |. 
although the inhomogeneous continent is sensitive to regularization effects, it is not obviously 
created by them. Indeed, there is no a priori reason to exclude the possibility of an inhomogeneous 
chiral symmetry breaking phase at high chemical potentials. For instance, inhomogeneous phases 
extending to arbitrarily high chemical potentials have been predicted for the Gross- Neveu model 



and its chiral counterpart [27|, for quarkyonic matter [221 ] and for QCD in the large- A'^c limit [18|. 
Here we do not want to enter this discussion, but simply take the inhomogeneous continent as 
a "model laboratory" to study the competition of one- and two-dimensional chiral crystals as a 
function of fi. 

Figure [5] shows the free energies associated with the one-dimensional solitonic solutions and the 
two-dimensional "egg carton" close to the phase transition marking the onset of the inhomoge- 
neous continent (the second-order nature of the transition from the restored phase is visible in the 
behaviour of both free energies). From this comparison, it is possible to see that in this region, 
two-dimensional solutions become favored over one-dimensional ones. Since with the current pa- 
rameter set the continent is not connected to the the inhomogeneous island, it is not clear at which 
point higher-dimensional structures become favored. In order to achieve a better understanding 
of the problem, in the following we therefore employ a slightly modified parameter set with a vac- 
uum constituent quark mass of 330 MeV. This does not modify any qualitative behaviour of the 
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FIG. 5. Free energies associated with the one-dimensional solitonic solutions (black solid line) and the 
two-dimensional "egg-carton" (blue dashed line) at the onset of the second inhomogeneous phase, occurring 
around /x = 500 MeV. Results are calculated for T = with a parameter set fitted to give a vacuum 
constituent quark mass of M^ac — 300 MeV, and are normalized with respect to the free energy of the 
restored phase. 
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FIG. 6. Free-energy diff'erences between inhomogeneous phases with different modulations, evaluated at 
T — Q and with a parameter set fitted to give a vacuum constituent quark mass of M^ac = 330 MeV. 



model but has the advantage that the inhomogeneous island merges with the continent, so that 
the comparison of the inhomogeneous phases can be performed on a continuous interval, without 
being interrupted by the restored phase. 

In Fig. [6] the differences between the free energies of three inhomogeneous phases with different 
modulations are displayed as functions of //. As we have seen before, at low chemical potentials 
the two-dimensional crystals are disfavored against the one-dimensional Jacobi elliptic function. 
Above /i ~ 450 MeV, however, the two-dimensional square lattice leads to a lower free energy. At 
fi w 600 MeV, also the hexagon surpasses the one-dimensional modulation and finally becomes the 
most favored shape at // ss 900 MeV. Thus, while being aware that the model cannot be trusted 
blindly in this high density region, it is nevertheless remarkable that we recover the same sequence 
of crystalline phases as described in Ref. 22l |. 
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V. DISCUSSION AND OUTLOOK 

In this paper we presented the results of our numerical study of two-dimensional chiral crystalline 
structures in the NJL model at zero temperature. At intermediate chemical potentials, in the region 
where the chiral phase transition would take place if the analysis was restricted to homogeneous 
condensates, we find that two-dimensional modulations are disfavored against one-dimensional 
ones, indicating that their greater kinetic energy cost is not sufficiently compensated by a larger 
gain in condensation energy. We also find that a hexagonal structure is even less favored than a 
square lattice in this regime. 

From these observations, it seems unlikely that a phase where the chiral condensate is modulated 
along three spatial dimensions could become thermodynamically favored, due to its even higher 
kinetic-energy cost. Moreover, since a GL analysis has revealed that also close to the Lifshitz point 



phases with higher-dimensional modulations are disfavored against one-dimensional ones [36(] , we do 
not expect that higher-dimensional structures appear at finite temperature, at least in mean-field 
approximation. A numerical analysis to confirm these expectations would of course be desirable. 

The situation gets successively reversed when we increase the chemical potential. At a certain 
value, we find that the structure of the ground state changes from having one-dimensional modula- 
tions to a two-dimensional square lattice and, at even higher ^u, to an hexagonal shape. Although 
the results must not be trusted blindly in this high-density regime, which lies at the edge of the 
expected range of validity of the model, we find it nevertheless noteworthy that we reproduce 
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It is also interes ting that 



qualitatively the behavior recently proposed for quarkyonic matter 

a similar sequence of phases was predicted for a 2D superconductor in a magnetic field [l2]. On 
the other hand, the chemical potentials where the two-dimensional structures appear in our model 
belong to the realm of color superconductivity. Therefore it would be important to include the 
effects of diquark pairing as well. 
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